/*
 * piece.cpp
 *
 *  Created on: 2012-12-5
 *      Author: shane
 */

#include "piece.h"
#include<cmath>
#include<typeinfo>
using namespace std;
namespace EME
{

simple_piece::simple_piece(const double& xmin_in, const double& xmax_in, const double& ymin_in, const double& ymax_in,
		const boundary& boundary_condition_in, const unsigned int& num_modes_in,
		const unsigned int& nx_in, const unsigned int& ny_in,
		pattern *const ptn_in) :
		piece(xmin_in, xmax_in, ymin_in, ymax_in, ptn_in, boundary_condition_in,
				num_modes_in,false), nx(nx_in), ny(ny_in)
{
	ex = vector<complex<double>*>(num_modes, NULL);
	ey = vector<complex<double>*>(num_modes, NULL);
	ez = vector<complex<double>*>(num_modes, NULL);
	hx = vector<complex<double>*>(num_modes, NULL);
	hy = vector<complex<double>*>(num_modes, NULL);
	hz = vector<complex<double>*>(num_modes, NULL);
	for (unsigned int i = 0; i < num_modes; ++i)
	{
		ex[i] = new complex<double> [nx * ny];
		ey[i] = new complex<double> [nx * ny];
		ez[i] = new complex<double> [nx * ny];
		hx[i] = new complex<double> [nx * ny];
		hy[i] = new complex<double> [nx * ny];
		hz[i] = new complex<double> [nx * ny];
	}
}
simple_piece::simple_piece(const simple_piece& B) :
		piece(B), nx(B.nx), ny(B.ny)
{
	ex = vector<complex<double>*>(num_modes, NULL);
	ex = vector<complex<double>*>(num_modes, NULL);
	ey = vector<complex<double>*>(num_modes, NULL);
	ez = vector<complex<double>*>(num_modes, NULL);
	hx = vector<complex<double>*>(num_modes, NULL);
	hy = vector<complex<double>*>(num_modes, NULL);
	hz = vector<complex<double>*>(num_modes, NULL);
	for (unsigned int i = 0; i < num_modes; ++i)
	{
		ex[i] = new complex<double> [nx * ny];
		ey[i] = new complex<double> [nx * ny];
		ex[i] = new complex<double> [nx * ny];
		hx[i] = new complex<double> [nx * ny];
		hy[i] = new complex<double> [nx * ny];
		hz[i] = new complex<double> [nx * ny];
		for (unsigned int j = 0; j < nx * ny; ++j)
		{
			ex[i][j] = B.ex[i][j];
			ey[i][j] = B.ey[i][j];
			ez[i][j] = B.ez[i][j];
			hx[i][j] = B.hx[i][j];
			hy[i][j] = B.hy[i][j];
			hz[i][j] = B.hz[i][j];
		}
	}
}
simple_piece::~simple_piece()
{
	for (unsigned int i = 0; i < num_modes; ++i)
	{
		delete[] ex[i];
		delete[] ey[i];
		delete[] ez[i];
		delete[] hx[i];
		delete[] hy[i];
		delete[] hz[i];
	}

}
double simple_piece::normalization(unsigned int n)
{
	if (!if_calculated)
		return 0.;
	double sum = 0;
	for (unsigned int i = 0; i < nx * ny; ++i)
	{
		sum += (ex[n][i] * conj(hy[n][i]) - ey[n][i] * conj(hx[n][i])).real()
				* 2.;
	}
	sum *= (xmax - xmin) / nx * (ymax - ymin) / ny;
	return sqrt(fabs(sum));
}
complex<double> simple_piece::projection(unsigned int na, signed char za,
		piece& B, unsigned int nb, signed char zb)
{
	if (typeid(B) == typeid(simple_piece))
		return projection(na, za, dynamic_cast<simple_piece&>(B), nb, zb);
	return 0.;
}
complex<double> simple_piece::projection(unsigned int na, signed char za,
		simple_piece& B, unsigned int nb, signed char zb)
{
	///point to point projection for the simplest situation
	if (nx == B.nx && ny == B.ny && xmin == B.xmin && xmax == B.xmax
			&& ymin == B.ymin && ymax == B.ymax)
	{
		complex<double> sum = 0.;
		for (unsigned int i = 0; i < nx * ny; ++i)
		{
			switch (zb)
			{
			case 1:
				sum += conj(ex[na][i]) * B.hy[nb][i]
						- conj(ey[na][i]) * B.hx[nb][i];
				break;
			case -1:
				sum -= conj(ex[na][i]) * B.hy[nb][i]
						- conj(ey[na][i]) * B.hx[nb][i];
				break;
			}
			switch (za)
			{
			case 1:
				sum += B.ex[nb][i] * conj(hy[na][i])
						- B.ey[nb][i] * conj(hx[na][i]);
				break;
			case -1:
				sum -= B.ex[nb][i] * conj(hy[na][i])
						- B.ey[nb][i] * conj(hx[na][i]);
				break;
			}
		}
		return sum / normalization(na) / B.normalization(nb);
	}
	///Otherwise the resolution and the position is not the same, we have more complicated case
	else
	{
		complex<double> sum = 0.;
		for (unsigned int i = 0; i < nx; ++i)
			for (unsigned int j = 0; j < ny; ++j)
			{
				rect_element a(xmin + (xmax - xmin) / nx * i,
						xmin + (xmax - xmin) / nx * (i + 1),
						ymin + (ymax - ymin) / ny * j,
						ymin + (ymax - ymin) / ny * (j + 1));
				double bdx = (B.xmax - B.xmin) / B.nx;
				unsigned int bi1, bi2;
				if (a.xmax <= B.xmin)
				{
					continue;
				}
				else if (a.xmin <= B.xmin)
				{
					bi1 = 0;
				}
				else
				{
					bi1 = static_cast<int>((a.xmin - B.xmin) / bdx);
				}
				if (a.xmin >= B.xmax)
				{
					continue;
				}
				else if (a.xmax >= B.xmax)
				{
					bi2 = nx - 1;
				}
				else
				{
					bi2 = nx - 1 - static_cast<int>((B.xmax - a.xmax) / bdx);
				}
				double bdy = (B.ymax - B.ymin) / B.ny;
				unsigned int bj1, bj2;
				if (a.ymax <= B.ymin)
				{
					continue;
				}
				else if (a.ymin <= B.ymin)
				{
					bj1 = 0;
				}
				else
				{
					bj1 = static_cast<int>((a.ymin - B.ymin) / bdy);
				}
				if (a.ymin >= B.ymax)
				{
					continue;
				}
				else if (a.ymax >= B.ymax)
				{
					bj2 = ny - 1;
				}
				else
				{
					bj2 = ny - 1 - static_cast<int>((B.ymax - a.ymax) / bdy);
				}
				for (unsigned int bi = bi1; bi <= bi2; ++bi)
					for (unsigned int bj = bj1; bj <= bj2; ++bj)
					{
						rect_element b(B.xmin + (B.xmax - B.xmin) / B.nx * bi,
								B.xmin + (B.xmax - B.xmin) / B.nx * (bi + 1),
								B.ymin + (B.ymax - B.ymin) / B.ny * bj,
								B.ymin + (B.ymax - B.ymin) / B.ny * (bj + 1));
						double area = a.overlap(b);
						complex<double> temp = 0.;
						unsigned int aa = i * ny + j;
						unsigned int bb = bi * B.ny + bj;
						switch (zb)
						{
						case 1:
							temp += conj(ex[na][aa]) * B.hy[nb][bb]
									- conj(ey[na][aa]) * B.hx[nb][bb];
							break;
						case -1:
							temp -= conj(ex[na][aa]) * B.hy[nb][bb]
									- conj(ey[na][aa]) * B.hx[nb][bb];
							break;
						}
						switch (za)
						{
						case 1:
							temp += B.ex[nb][bb] * conj(hy[na][aa])
									- B.ey[nb][bb] * conj(hx[na][aa]);
							break;
						case -1:
							temp -= B.ex[nb][bb] * conj(hy[na][aa])
									- B.ey[nb][bb] * conj(hx[na][aa]);
							break;
						}
						temp *= area;
						sum += temp;
					}
			}
		return sum / normalization(na) / B.normalization(nb);
	}
	return 0.;
}
}
